/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TAEikonalPhase
  \brief compute eikonal phase \chi=\chi_N + \chi_C, S=exp(i*\chi), where \chi_N
  represents nuclear phase, and \chi_C Coulomb phase.
  This is a class.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/07/25
  \date Last modified: 2020/09/05 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#ifndef _TAEikonalPhase_h_
#define _TAEikonalPhase_h_

#include <complex>
#include <string>

using std::string;
typedef std::complex<double> cdouble;

class TAEikonalPhase{
public:
  // alphaP, alphaT: proj(targ) nucleon size factor: rho_n(r)=exp(-r^2/alphap^2)
  TAEikonalPhase(int zP, int aP, int zT, int aT, double ek, double alphaP, double alphaT);
  virtual ~TAEikonalPhase();

  /// \retval eikonal nucclear phase = 1/KNN*\int_0^\infty{q*rhoP*rhoT*fNN*J0(qb)}
  /// \param b: impact parameter
  virtual cdouble GetPhaseN(double b);
  virtual double GetPhaseNI(double b); ///< \return the imag part of the phase
  virtual double GetPhaseNR(double b); ///< \return the real part of the phase
  /// \retval eikonal Coulumb phase = 2\eta*ln(kb)
  /// \param b: impact parameter
  virtual double GetPhaseC(double b) const;
  /// \retval \chi_C+\chi_N
  virtual cdouble GetPhase(double b);
  double GetAlphaNN() const{ return fAlphaNN; }
  /// Fourier transform of nucleus density on the transverse plane
  /// alphap is the width (sigma) of proton (nucleon) Gaussian density
  void ReadDensityFile(const string &file, int &n, double *&r, double *&rho, int a);
  void SetTargetDensity(const string &file);
  void SetProjectileDensity(const string &file);
  /// set optical potentials
  void SetOpticalPotential(const string &opt);
  double PM() const{ return fPM; }
  double TM() const{ return fTM; }

  static const double QMAX; ///< the maximum momentum transfer (in fm^-1)

  friend class TASMatrix;

protected:
  double fZP, fAP, fZT, fAT; ///< projectile and target (N, Z)
  double fEk; ///< Elab in MeV/nucleon
  double fMu; ///< the reduced mass in MeV (mP*mT)/(mP+mT)
  /// Array R in those densities are supposed to be equidistant
  int fNRP, fNRT; ///< length of fRP and fRT
  double *fRP, *fRhoP; ///< density of the projectile (fRP, fRhoP) in r (fm)
  double *fRT, *fRhoT; ///< density of the target (fRT, fRhoT) in r (fm)
  double fAlphaP, fAlphaT; ///< ap: exp(-r^2/ap^2) for projectile and target
  double *fQ, *fRhoQ; ///< = q*rhoP(q)*rhoT(q)

  double fAlphaNN, fBetaNN, fSigmaNN;
  double fBeta; ///< the relative speed: P w.r.t. T
  double fPM, fTM; ///< mass in MeV of the projectile and the target
  bool fIsPauli; ///< true: consider Pauli blocking; false: otherwise

  bool fIsOpt; ///< if using an external optical potential, instead of t-rho-rho
  int fOptN; // length of the input optical potential arrays
  double *fOptR, *fV, *fW; ///< the optical potential: (fV+i*fW)(fOptR)
};

#endif
